Innovative GenExpA software for selecting suitable reference genes for reliable normalization of gene expression in melanoma

The algorithms commonly used to select the best stable reference gene in RT-qPCR data analysis have their limitations. We showed that simple selection of the reference gene or pair of genes with the lowest stability value from the pool of potential reference genes—a commonly used approach—is not sufficient to accurately and reliably normalize the target gene transcript and can lead to biologically incorrect conclusions. For reliable assessment of changes in a target gene expression level, we propose our innovative GenExpA software, which works in a manner independent of the experimental model and the normalizer used. GenExpA software selects the best reference by combining the NormFinder algorithm with progressive removal of the least stable gene from the candidate genes in a given experimental model and in the set of daughter models assigned to it. The reliability of references is validated based on the consistency of the statistical analyses of normalized target gene expression levels through all models, described by the coherence score (CS). The use of the CS value imparts a new quality to qPCR analysis because it clarifies how low the stability value of reference must be in order for biologically correct conclusions to be drawn. We tested our method on qPCR data for the B4GALT genes family in melanoma, which is characterized by a high mutation rate, and in melanocytes. GenExpA is available at https://github.com/DorotaHojaLukowicz/GenExpA or https://www.sciencemarket.pl/baza-programow-open-source#oferty.

Achieving reliable results for normalization of the transcript level of a target gene requires an internal reference gene or pair of genes (housekeeping genes, HKGs) with stable expression between the analyzed samples and under different experimental conditions. Moreover, it is believed that the transcript levels of the reference and target genes should be expressed at roughly the same level 1 . However, the transcript levels of many of the HKGs typically used as references may change significantly under physiological or pathological conditions as well as across tissue types and experimental conditions [1][2][3][4][5][6][7][8] . Therefore the selection and validation of the normalizer are the most critical stages in qPCR data analysis. Typically, one or more of five algorithms (geNorm 3 , NormFinder 9 , BestKeeper 10 , the comparative ΔCt method 11 and RefFinder; http:// 150. 216. 56. 64/ refer enceg ene. php) are used to select the best stable single or combination of reference genes from a panel of candidate genes, although there are also other proposed algorithms 12 . The normalizer with the lowest stability value is considered the best and is then used for calculation of target gene expression in a group of samples. We have recently shown that this commonly accepted approach seems not enough for accurate and reliable target gene transcript normalization and that it may lead to biologically incorrect conclusions 13 . We demonstrated that beyond the parent model of samples of interest (experimental model), it is necessary to construct daughter models (auxiliary models built from the same samples as the experimental model but with fewer samples-combinations of samples of interest without repetition), followed by selection and validation of an appropriate normalizer for each model. Here we have shown that validation of the normalizer is done by determining the coherence score for target gene expression analyses performed for all given models (experimental and daughter models). If the statistical analysis produces inconsistent results for the expression level of a target gene based on a comparison of individual sets of two samples through all models, we have to search for new references with improved stability values via progressive removal of the least stable candidate reference gene in each sample, followed by re-selection and www.nature.com/scientificreports/ re-validation of new normalizers each time 13 . However, this proposed approach of RT-qPCR data analysis has been a very time-consuming task, requiring the use of several specialist software packages. This problem led us to automate these calculations by developing the GenExpA tool (Gene Expression Analyzer). It is a comprehensive tool based on a previously described workflow for quantified as well as raw qPCR data 13 . Based on the selected reference gene or pair of genes, GenExpA calculates the relative target gene expression level, performs statistical analyses and generates results in the form of graphs and tables. The innovative aspect of GenExpA is that it estimates the coherence of results for each target gene in all designed models. It determines the robustness of normalization and works out the coherence score for the individual target gene. The advantage of this strategy is its ability to determine relative target gene expression based on the consistency of the results in all tested samples, regardless of the experimental model and the reference gene or pair of genes used. Here, based on qPCR data for B4GalT1-B4GalT7 transcripts in melanocytes and melanoma cells, we show how to perform the analysis step by step, and how to solve problems using this new methodological approach.  In the setting mode (not shown in this paper; see the GenExpA manual), the user separates 'n' potential reference genes from a list of all genes. After clicking on 'Generate combinations' , the GenExpA tool automatically creates a set of possible 'z' models, i.e., 'z' models composed of samples of interest, which are combinations without repetitions (order does not matter) of two or more samples from the pool of 'x' samples. The 'z' models consist of the experimental model and its 'z − 1' daughter models (auxiliary models); for example, if the experimental model consists of 5 samples, the number of daughter models is 25 (10 models of 2 samples, 10 models of 3 samples, and 5 models of 4 samples). The user unchecks the statistical model, wherein the pairwise t-test with Holm adjustment is recommended as a statistical model for multiple pairs of observations in the case of a normal distribution (to determine the significance of a difference in gene expression between two or more samples). Alternatively, the non-parametric Mann-Whitney test or the ANOVA Kruskal-Wallis test with subsequent post hoc Dunn's test are recommended for models composed of two unmatched samples or three or more unmatched samples, respectively, in the case of a non-normal distribution 23 . In the setting mode, the user determines the P value indicating the level of statistical significance (default P < 0.05; not shown in this paper; see the GenExpA manual). The user starts the automated data analysis by clicking on 'Run calculations'; it causes the implemented NormFinder algorithm to select, in each of the analyzed models, the best references between 'n' potential reference genes (purple line) according to the stability values of their expression in the experimental model and its auxiliary models. The best stability value is referred to the minimal combined inter-and intrasample expression variation of the gene 9 . NormFinder works on Ct values converted to quantified values, which were calculated based on the calibration curve previously prepared for all candidate reference genes. Then the GenExpA tool normalizes the expression of target genes (RQ value) in each model in relation to selected references, conducts a statistical analysis and calculates the coherence score (CS) separately for each target gene [details of how GenExpA determines CS are shown in (b)]. If the CS value is equal to 1, the target gene expression analysis is completed. In the case of inconsistent results (CS < 1), the user can change the setting parameters by choosing 1 in the 'Remove repetitions' box and by marking the option 'Select best remove for models' (not shown in this paper; see the GenExpA manual). These new settings cause the least stable HKG in each model to be removed from the pool of 'n' candidate reference genes. The user starts the improved analysis by clicking on 'Run calculations' . The GenExpA tool chooses new references based on a reduced pool of HKGs to 'n − 1' in each model and conducts a reanalysis (statistic and CS calculation) of the normalized target gene(s) (red line). This approach can be repeated serially (number in 'Remove repetitions' box should be increased by 1 each time) when the CS value is below 1 and the pool of reference genes is not less than three (from the green line to the navy blue line). It is important to note that after marking the option ‛Select best remove for models' , GenExpA selects the reference for a given model from the level at which the stability value was the lowest.

Results
A flow chart of target gene expression and statistical analysis of qPCR data using the GenExpA tool is shown in Fig. 1. At the beginning, we uploaded raw qPCR data for candidate reference (HPRT1, PGK1, RPS23, SNRPA; set 1) and target (B4GALT1-B4GALT7) genes (Supplementary Table S1) as well as quantified qPCR data for candidate reference genes (Supplementary Table S2) obtained for five human cell lines (melanocytes and melanoma cells from different stages of oncogenic progression) comprising the experimental model. It is important to note that here "cell line" means "sample". Since all PCR reactions were performed in three technical repeats for every Table 1. Reference genes and their stability values in a panel of the experimental model and daughter models selected by NormFinder software on the base of 4 (part A) and 5 (part B) potential reference genes before (0) and after removing one (1) or two (2) the least stable genes from a given set of 4 and 5 reference genes. The differences in the choice of the normalizer or obtained stability values are written in bold. www.nature.com/scientificreports/ creates an opportunity to select the reference gene with a lower stability value in the experimental model as well as auxiliary models, and therefore gives a more robust analysis of target gene expression. For example, the use of a set of five candidate reference genes (PGK1, RPS23, SNRPA, HPRT1, GUSB) with sequential removal of the two least stable potential reference genes gave the lowest stability value of the selected reference gene for the experimental model (0.215) and auxiliary models (0.010-0.239). On the other hand, we noticed that box-plots 13 and 91 show exactly the same results as, respectively, box-plots 2, 5, 6, 10, 12 and 80, 83, 84, 90 (compare the medians of relative quantity (RQ) and statistical significance between samples in Supplementary Fig. S19). It means that in the case of B4GALT1 and B4GALT7 gene expression analysis, the use of four HKGs from sets 2 or 5 with rejection of the most unstable HKG, or from set 3, is sufficient. The same is true for genes B4GALT2 and B4GALT3 and set 3 (compare box-plots 26 and 39 with 18, 19 and 31, 32, respectively, in Supplementary  Fig. S19). Similar relationships occur for gene B4GALT5 and sets 3 and 1 after rejection of the most unstable HKG (compare box-plots 78 with 57, 58 and 54 in Supplementary Fig. S19). Therefore, adding a fifth housekeeping gene is mandatory to confirm the robustness of an analysis conducted with the above-mentioned sets of four candidate reference genes. Interestingly, the lower the HKG stability value, the better the CS value, that is, the more robust the analysis, potentially much more correct biologically.
On the other hand, selection of the normalizer based on five potential reference genes and with 'Remove repetition' set at 2 led to an inconsistent analysis for the relative expression of gene B4GALT6 (box-plot 78 in Supplementary Fig. S19), although the same box-plots were generated with 'Remove repetition' set at 1 as well as in the case of set 3 with 'Remove repetition' set at 0 or 1 (compare box-plots 70, 71, 77 and 78 in Supplementary  Fig. S19), and the reference stability values for the experimental model remained the same (0.215). It means that if a more robust analysis is needed, to improve its consistency another HKG/s should be introduced into the pool of candidate reference genes. Another approach may be to complete the analysis with 'Remove repetition' set at 1 (box-plot 77 in Supplementary Fig. S19) with indication of the CS values and stability parameters for the reference genes selected in the experimental and auxiliary models.
We want to emphasize that simple application of the NormFinder algorithm ('Remove repetition' 0) may lead to a biological misinterpretation of the obtained results, as it strongly depends on the choice of a set of potential reference genes. For example, sets 1 and 5 showed significantly higher relative expression of gene B4GALT6 in WM793 cells than in WM266-4 cells. In turn, sets 2, 4 and the set of five HKGs showed, in these two cell lines, an inverse relationship between the B4GALT6 gene expression that was also statistically significant (compare box-plots 66 and 74 with box-plots 68, 72 and 76 in Supplementary Fig. S19). Moreover, the relative expression levels of B4GALT5, B4GALT6 and B4GALT 7 were significantly lower in cell line WM35 than in cell line WM266-4, as shown in sets 2 and 5 of candidate reference genes (compare box-plots 55, 63 and 68, 76 and 81, 89 with 59 and 72 and 85 in Supplementary Fig. S19), or significantly higher as shown in set 4 (box-plot 72 in Supplementary Fig. S19).
Confusing results are also observed for gene B4GALT3: in sets 1 (box-plot 27 in Supplementary Fig. S19), 2 (box-plot 29 in Supplementary Fig. S19), 3 (box-plot 31 in Supplementary Fig. S19), 5 (box-plot 35 in Supplementary Fig. S19) and in set of five candidate reference genes (box-plot 37 in Supplementary Fig. S19) its relative expression level was significantly lower in cell line WM35 than in cell line WM793, but in set 4 the relationship was the reverse (box-plot 33 in Supplementary Fig. S19). Similar for the gene B4GALT6: in sets 1 (box-plot 66 in Supplementary Fig. S19) and 5 (box-plot 74 in Supplementary Fig. S19) its relative expression level was significantly lower in cell line WM35 than in cell line WM793, but in set 4 the relationship was the reverse (boxplot 72 in Supplementary Fig. S19). In addition, a larger pool of candidate reference genes does not always lead to selecting the better normalizer in simple application of the NormFinder algorithm (compare the reference stability values for the experimental model and auxiliary models in set 3 and the set of five candidate reference genes; 'Remove repetition' 0). In contrast, our method based on gradual removal of the gene with the highest variability of expression leads to selection of a normalizer with a lower stability value (compare the reference stability values for the experimental model and auxiliary models in sets 1-5 with 'Remove repetition' 1 and the set of five candidate reference genes with 'Remove repetition' 2). Figure 2 shows the results of GenExpA analysis for the expression levels of target genes B4GALT1-B4GALT7 in the analyzed experimental model. We obtained consistent and complete analyses for all target genes under the given stability values for the experimental and auxiliary models. These presented results demonstrate that in melanoma cells, in particular metastatic melanoma cell line WM266-4, the expression of all B4GALT genes decreases relative to that of melanocytes. An exception is melanoma cell line WM793, in which the expression levels of genes B4GALT2 and B4GALT4 show no significant differences as compared to those for melanocytes.

Discussion
Obtaining reliable results from the PCR reaction is the outcome of many factors related to handling of the material, beginning with the step of RNA/DNA isolation and ending in analysis of the results. That is why the MIQE guidelines were proposed for transparent reporting of experimental details in all publication prepared with the use of qPCR technique 14 . The reaction needs to be standardized in molecular analyses of all kinds of biological materials [15][16][17][18] . Over the past decade the MIQE guidelines have not come into common use. Dijkstra et al. 15 made a critical analysis of qPCR result normalization in 179 publications regarding colorectal cancer. They showed that only 3% of these studies applied qPCR methodology based on the use and validation of multiple reference genes. Several statistical are used to select the best stable single gene or combination of reference genes from a panel of candidate genes. It should be stressed that each of these algorithms has its limitations which affect the ranking of candidate reference genes and finally the selection of the best normalizer 19,20 . Recent studies have proposed the use of at least three different methods, but the problem of how to interpret conflicting results obtained from the use of different statistical methods still puzzles researchers. They try to integrate a few algorithms to average the www.nature.com/scientificreports/  Table 1, part B). The statistical analyses used the pairwise t-test with Holm adjustment. Red line represents statistical significance, P < 0.05. www.nature.com/scientificreports/ stability ranks; for example, the RefFinder algorithm calculates the geometric mean of ranking values obtained by four other algorithms 19 . But these statistical methods differ from each other; averaging the ranks can lead to a suboptimal assessment of stable reference genes 20 . It is commonly accepted that the lower the stability value of the reference, the greater the certainty of correct determination of the relative expression of the target gene. Yet it has not been clarified how low the stability value must be. In this paper we showed that the procedure for choosing a suitable reference involves two steps: applying an algorithm for best reference gene selection (here, NormFinder), and checking whether the stability level of the selected reference is low enough to lead to a correct biological interpretation of relative target gene expression (here, determination of CS for the analysis). Through the application of these two steps, our method validates the chosen reference based on determination of the coherence score for the normalized target gene. In our method we used the model-dependent NormFinder algorithm, and by simulating the removal of the least stable gene from a set of candidate genes in the experimental as well as daughter models we were able to obtain more stable normalizers for each model. Previously, this approach helped us to exclude incorrect normalizations and biological misinterpretations concerning the alteration of ERAP1 gene expression during melanoma progression 13 . We chose the NormFinder algorithm because it calculates the stability of a reference gene or pair of genes by analysing the variation of its expression within and between samples. However, genes with high overall variation influence the stability ranking of candidate reference genes. Herein we have demonstrated that simple application of NormFinder ('Remove repetition' 0) may lead to biological misinterpretation of the results, because it depends on the HKG set used to select the reference. Our method eliminates the influence of the most variable gene or genes on the stability rankings of all other candidate reference genes. In addition, breaking down the experimental model into daughter models (auxiliary models) allows us to check whether a uniform trend of target gene expression among particular cell lines is maintained. This trend decides the consistency of the analysis and is described by coherence score CS = 1.
Obtaining full consistency of the analysis is tantamount to its completion. As we suggest here, based on our B4GALT gene expression analysis, the lower the stability value, the better the coherence of the obtained results, which, we propose, is related to their biological correctness. If, despite the application of our strategy, a uniform trend of target gene expression among particular cell lines has not been achieved, introducing a successive HKG/s, along with repeated rejection of the weakest gene, can lead to selection of the better reference and, in turn, biologically correct conclusions. It is important to note that once selected, a reference for a given model will not necessarily be selected once again in an experiment repeated for this model after some time, even if the pool of candidate reference genes used is exactly the same. This is because the gene expression pattern at the time point of sample collection depends on the current rate of cell mass growth and on microenvironmental conditions 21,22 .

Conclusions
Here we have presented an expanded version of our already published method for RT-qPCR data analysis implemented in the GeneExpA tool-software for reference gene validation and automatic calculation of target gene expression. GenExpA will assign a CS value ranging from 0 to 1 to each analysis; a fully coherent analysis is described by a CS value of 1. We suggest starting the analysis with a set of at least four potential HKGs. If the results are not satisfactory-if the CS value is below 1-introducing another candidate HKG/s to the set of candidate reference genes may improve them. However, even if CS is equal to 1 for a given target gene, enlarging the pool of candidate reference genes can lead to an analysis that is more robust in terms of between-sample statistical significance in the experimental model. Interestingly, our workflow, including full target gene analysis in order to validate normalizer, shows that the lower stability value of the selected reference gene or pair of genes is related to the biological correctness of the results. Moreover, our workflow is the first one that allows the user to define the required stability value needed to draw biologically correct conclusions. For clarity of the biological conclusions drawn, the description of the results should give stability values assigned to the experimental model www.nature.com/scientificreports/ and auxiliary models; this may help other researchers to compare their results and understand any differences between their outcomes and yours. GenExpA software can carry out an analysis of many genes independently at the same time. RT-qPCR data analysis. The GenExpA tool was used for fully automated qPCR data analysis. The gene or combination of two HKGs with the lowest stability value were selected as the best reference by combining the model-based variance estimation with progressive removal of the least stable reference gene. For statistical analysis, the pairwise t-test with Holm adjustment was used. P values less than 0.05 were taken to indicate statistical significance (P < 0.05).